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Abstract 

The  dynamic  thermomechanical  response  of  a  tungsten  heavy  alloy  is  investigated  via  modeling  and  numerical  sim¬ 
ulation.  The  material  of  study  consists  of  relatively  stiff  pure  tungsten  grains  embedded  within  a  more  ductile  matrix 
alloy  comprised  of  tungsten,  nickel,  and  iron.  Constitutive  models  implemented  for  each  phase  account  for  finite  defor¬ 
mation,  heat  conduction,  plastic  anisotropy,  strain-rate  dependence  of  flow  stress,  thermal  softening,  and  thermoelastic 
coupling.  The  potentially  nonlinear  volumetric  response  in  tungsten  at  large  pressures  is  addressed  by  a  pressure-depen¬ 
dent  effective  bulk  modulus.  Our  framework  also  provides  a  quantitative  prediction  of  the  total  dislocation  density, 
associated  with  cumulative  strain  hardening  in  each  phase,  and  enables  calculation  of  the  fraction  of  plastic  dissipation 
converted  into  heat  energy.  Cohesive  failure  models  are  employed  to  represent  intergranular  fracture  at  grain  and  phase 
boundaries.  Dynamic  finite  element  simulations  illustrate  the  response  of  realistic  volume  elements  of  the  poly  crystal¬ 
line  micro  structure  subjected  to  compressive  impact  loadings,  ultimately  resulting  in  spallation  of  the  material.  The  rel¬ 
ative  effects  of  mixed-mode  interfacial  failure  criteria,  thermally-dependent  fracture  strengths,  and  grain  shapes  and 
orientations  upon  spall  behavior  are  weighed,  with  interfacial  properties  exerting  a  somewhat  larger  influence  on  the 
average  pressure  supported  by  the  volume  element  than  grain  shapes  and  initial  lattice  orientations  within  the  bulk 
material.  Spatially  resolved  profiles  of  particle  velocities  at  the  free  surfaces  of  the  volume  elements  indicate  the  degree 
to  which  the  incident  and  reflected  stress  waves  are  altered  by  the  heterogeneous  micro  structure. 
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1.  Introduction 

Because  tungsten  heavy  alloys  (WHAs)  typically  exhibit  a  relatively  large  mass  density,  high  melting 
point,  and  high  failure  strength  at  elevated  rates  of  loading,  these  materials  are  appealing  for  implementa¬ 
tion  in  long-rod  kinetic  energy  penetration  applications  (Magness,  1994;  Cai  et  al.,  1995).  The  particular 
material  investigated  here  consists  of  pure  tungsten  grains  (BCC)  embedded  in  a  polycrystalline  matrix 
phase  (FCC)  comprised  of  nickel  (50  wt.%),  iron  (25  wt.%),  and  tungsten  (25  wt.%).  Fabrication  of  the 
composite  material  is  conducted  by  isostatic  pressing  and  sintering  of  a  mixture  of  W,  Ni,  and  Fe  powders, 
followed  by  an  annealing  step  to  remove  absorbed  hydrogen  and  then  possible  quenching,  swaging,  and/or 
pre-twisting  to  alter  the  dynamic  mechanical  properties  prior  to  deformation  testing  at  high  strain  rates 
(Weerasooriya,  1998;  Wei  et  ah,  2000).  Typical  grain  sizes  are  10-30  pm  for  the  W  crystals  and  200- 
500  pm  for  the  matrix  phase  (Zhou,  1993),  meaning  that  multiple  W  crystals  are  often  embedded  within 
each  single  “grain”  of  the  matrix  phase. 

Macroscopic  constitutive  models  (cf.  Stevens  and  Batra,  1998;  Wei  et  al.,  2001)  previously  used  to  rep¬ 
resent  the  two-phase  system  in  numerical  simulations  lack  a  rigorous  description  of  the  kinematics  of  aniso¬ 
tropic  plasticity  (e.g.  plastic  spin),  the  role  of  lattice  orientations  of  constituent  grains,  and  the  effects  of 
crystal  morphology  on  failure  processes.  Microstructural  characteristics  such  as  grain  shape,  spatial 
arrangement  of  grains,  and  local  crystallographic  orientation  are  known  to  influence  the  material’s  thermo¬ 
mechanical  response.  Furthermore,  experimental  (Bruchey  et  al.,  1991,  1992)  and  numerical  (Schoenfeld 
and  Benson,  1997)  studies  have  suggested  a  possible  correlation  between  the  lattice  orientation  of  W  single 
crystals  and  the  performance  of  such  crystals  as  kinetic  energy  penetrators.  Previous  computational  work 
by  Zhou  and  co-workers  (Zhou  et  al.,  1994;  Zhou,  1998)  explained  the  influences  of  grain  morphology,  heat 
conduction,  strain-  and  strain-rate  hardening,  and  thermal  softening  on  the  elastoplastic  deformation  and 
shear  localization  of  tungsten  alloys  at  high  rates  of  shear  and  combined  pressure-shear  loading.  While 
grain  and  matrix  phases  were  resolved  discretely  in  these  calculations,  each  constituent  was  treated  as  an 
isotropically-hardening,  hypothermoelastic-viscoplastic  material,  thus  neglecting  possibly  significant  effects 
of  lattice  rotation  and  anisotropic  strain  hardening  associated  with  the  evolution  of  crystallographic  tex¬ 
ture.  In  contrast,  in  the  present  investigation  we  employ  a  crystal  plasticity  framework  (cf.  Asaro,  1983) 
to  examine  the  effects  of  crystallographic  orientation  and  plastic  anisotropy. 

When  subjected  to  tensile  states  of  stress,  tungsten  alloys  are  known  to  fracture  readily  along  grain  and 
phase  boundaries  (Zamora  et  al.,  1992;  Woodward  and  O’Donnell,  2000;  Weerasooriya,  2003).  Wee¬ 
rasooriya  (2003)  deduced  that  failure  of  WHA  specimens  subjected  to  high  rates  of  tensile  deformation  of¬ 
ten  initiates  via  fracture  at  W-W  interfaces,  and  less  often  at  interfaces  between  W  grains  and  the  matrix 
phase.  The  former  interfaces  are  often  thought  to  be  weaker  than  the  latter,  based  upon  visual  examination 
of  recovered  specimens  that  demonstrate  how  fracture  surfaces  are  prone  to  initiate  at  W-W  boundaries 
and  then  either  propagate  along  W-matrix  interfaces,  or  less  frequently,  along  cleavage  planes  within  grains 
of  either  phase  (Zamora  et  al.,  1992;  Weerasooriya  et  al.,  1994;  O’Donnell  and  Woodward,  2000;  Wood¬ 
ward  and  O’Donnell,  2000;  Weerasooriya,  2003).  O’Donnell  and  Woodward  (2000)  noted  that  at  higher 
ambient  temperatures  (>373  K),  the  tendency  for  cleavage  fracture  within  grains  tends  to  increase  at  the 
expense  of  interfacial  fracture.  Zamora  et  al.  (1992)  observed  that  microscopic  fracture  mechanisms  and 
macroscopic  fracture  toughness  measurements  in  such  alloys  vary  widely  with  chemical  composition,  grain 
size,  and  grain  continuity.  In  our  constitutive  framework  and  associated  numerical  scheme,  we  address 
interfacial  fracture  through  cohesive  zone  models  (Barenblatt,  1959;  Camacho  and  Ortiz,  1996).  Thus  we 
are  able  to  capture  a  failure  mode  that  may  dominate  the  response  of  the  WHA  system  under  tensile  stress 
states,  in  contrast  to  prior  computational  studies  reported  in  the  literature  that  did  not  treat  fracture 
explicitly. 

In  the  present  work  we  investigate  tensile  fractures  emerging  in  the  setting  of  a  numerical  wave  propa¬ 
gation  simulation  (cf.  Espinosa  and  Zavattieri,  2003a, b).  A  finite  sample  of  the  material  is  subjected  to  a 
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compressive  pressure  pulse  which  subsequently  traverses  the  sample  and  then  reflects  off  of  its  free  surface 
as  a  tensile  wave,  leading  to  spall  failure  within  the  sample.  Dandekar  (1976)  and  Asay  et  al.  (1980)  studied 
the  shock  behavior  of  pure  polycrystalline  tungsten.  The  high  pressure  responses,  including  spallation,  of 
various  W  alloys  were  probed  experimentally  via  plate  impact  testing  by  Dandekar  and  Weisgerber  (1999) 
and  through  explosive  loading  by  Baoping  et  al.  (1994).  The  former  researchers  (Dandekar  and  Weisgerber, 
1999)  found  the  spall  stress  of  WHAs  to  be  nearly  independent  of  impact  stress  over  a  range  of  5.9-24  GPa, 
attaining  values  of  1.9  ±  0.4  GPa,  whereas  the  latter  (Baoping  et  al.,  1994)  reported  an  increase  in  spall 
strength  of  WHAs  from  0.54  GPa  to  3  GPa  when  shocked  to  applied  stress  levels  between  21  and 
30  GPa,  followed  by  a  decrease  in  strength  to  1.56  GPa  at  an  applied  pressure  of  32  GPa.  Baoping  et  al. 
(1994)  noted  that  atomic  diffusion  or  phase  changes  may  affect  the  response  of  the  material  at  extremely 
high  pressures,  as  may  a  transition  in  dominant  failure  mechanism  from  intragranular  to  intergranular  frac¬ 
tures,  or  vice-versa. 

The  present  study  employs  single  crystalline  plasticity  models  for  the  bulk  response  of  each  phase,  as 
detailed  in  Section  2.  Constitutive  models  for  both  materials  are  embedded  within  a  general  thermodynamic 
framework  applicable  for  describing  any  thermoelastic-viscoplastically-deforming  metallic  crystal.  The 
description  accounts  for  finite  deformation  and  volume  change,  strain-rate  dependence,  heat  conduction, 
thermal  expansion,  thermal  softening,  and  thermoelastic  coupling.  Furthermore,  our  framework  enables 
calculation  of  the  percentage  of  plastic  work  stored  in  the  material  as  residual  elastic  energy  associated  with 
lattice  defects  (i.e.,  dislocations).  As  discussed  in  Section  3,  cohesive  zone  models  accounting  for  failure 
strengths  depending  upon  mode  mixity  and  temperature  represent  the  constitutive  behavior  in  damaged  re¬ 
gions  at  intergranular  interfaces.  Finite  element  meshes  of  equiaxed  grains  and  elongated  grains  are  con¬ 
structed,  corresponding  to  transverse  and  longitudinal  sections,  respectively,  of  samples  of  WHA 
material  cut  from  an  actual  kinetic  energy  penetrator.  Experimentally-measured  crystallographic  textures 
and  randomly-generated  textures  are  also  varied  as  initial  conditions  in  the  simulations.  Poly  crystals  are 
subjected  to  high-rate  uniaxial  compression  leading  to  spall  fracture,  with  results  of  these  simulations  de¬ 
scribed  in  Section  4. 


2.  Crystal  plasticity  formulation 

2.1.  Kinematics 


Let  x  =  x(X,  t )  represent  a  smooth  time-dependent  motion,  with  x  and  X  spatial  and  reference  coordi¬ 
nates,  respectively.  The  local  deformation  gradient  f  is  decomposed  as 


f  = 


0X 

ax 


rr  °fp, 


(i) 


where  F,  f0,  and  F  represent,  respectively,  the  kinematics  of  elasticity  and  rigid-body  rotation,  stress-free 
thermal  expansion  or  contraction,  and  the  cumulative  contribution  of  moving  crystal  defects.  In  Eq.  (1) 
and  henceforward,  juxtaposition  implies  summation  over  one  set  of  adjacent  indices,  i.e.,  (AB)  ^  =  AljBJk 
for  any  second-rank  tensors  A  and  B.  During  arbitrarily  heterogeneous  plastic  flow,  none  of  F,  f6,  and 
F  is  a  compatible  (i.e.,  integrable)  deformation  gradient.  As  shown  in  Fig.  1,  the  elastic  and  thermal  terms 
dictate  the  deformation  of  the  slip  direction  contravariant  vectors  s(a)  and  slip  plane  normal  covariant  vec¬ 
tors  m(a): 


s(a)  =  rf0s(«)? 


m(a)  = 


\re-\ 


(2) 


with  the  superposed  1”  denoting  matrix  inversion.  The  spatial  velocity  gradient  1  is  written  as 
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Fig.  1.  Deformations,  configurations,  and  slip  vectors. 


0X 

0X  : 


ff“ 


1  _|_  ^e£9^0—  1  1  _|_  ^e^0^p^p—\^9 


-1  ire— 1 


=F 


(3) 


where  the  superposed  dot  implies  a  material  time  derivative.  As  we  are  concerned  here  with  cubic  lattices, 
thermal  deformation  is  assumed  isotropic  (Lee  et  ah,  1997) 

1°  =  fV”1  =  «T01,  (4) 

where  6  is  the  temperature  change  measured  from  the  reference  state  and  aT(0)  is  the  temperature-depen¬ 
dent  thermal  expansion  coefficient  giving  the  change  in  length  per  unit  current  length  per  unit  increment  in 
6.  The  unit  tensor  is  denoted  by  1.  The  plastic  velocity  gradient  in  the  intermediate  configuration  b  of  Fig.  1 
is  defined  as  in  classical  crystal  plasticity  theory  (Asaro,  1983): 

F  =  fPfP- 1  =  7(a)sf  <g>  ,  (5) 

a=l 

with  y ^  the  plastic  shearing  rate  on  slip  system  a,  n  the  number  of  potentially  active  slip  systems,  and  (g)  the 
tensor  product.  As  a  result  of  the  orthogonality  of  and  m^,  the  plastic  deformation  is  volume-preserv¬ 
ing,  i.e.  trP  =  0,  where  the  trace  operation  for  a  second  rank  tensor  A  is  denoted  by  tr(A)  =  Akk. 


2.2.  Balance  laws 


Localized  balance  laws  are  specified  in  the  current  configuration  as  follows: 

p  +  ptr(l)  =  0,  div<7  +  pb  =  px,  a  =  <rT,  pe  +  divq  —  tr(o-gl)  =  pr,  (6) 

with  p,  g ,  b,  e ,  q,  and  r  denoting  the  current  mass  density,  contravariant  Cauchy  stress  tensor,  body  force 
vector  per  unit  mass,  internal  energy  per  unit  mass,  heat  flux  vector  per  unit  current  area,  and  energy  source 
per  unit  mass,  respectively.  We  use  ‘div’  to  represent  divergence  with  respect  to  spatial  coordinates,  for 
example  divq  =  q“a  =  daqa  +  yaabqb  with  the  usual  Christoffel  symbols  in  the  spatial  coordinate  system  sat¬ 
isfying  2 gadydhc  —  dcgba  +  dbgca  ~  0 agcb ?  where  the  compact  notation  da  =  d/dxa  and  the  spatial  metric  satis¬ 
fies  gab  =  0aX'0^x.  The  local  entropy  inequality  is  written  as  follows,  with  r\  the  time  rate  of  entropy 
production  per  unit  mass 

q 


pr 
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Next  introduced  is  the  Helmholtz  free  energy  per  unit  mass,  \[r. 


\j/  =  e  —  Orj. 


(8) 


Subsequently,  upon  substitution  of  (6)  and  (8)  into  (7),  the  entropy  inequality  becomes 

tr(ffgl)  -  q  S*  p{\ \  +  Ori),  (9) 

with  Vx0  the  contravariant  temperature  gradient  in  the  spatial  frame  and  •  a  scalar  product  operation  for 
vectors,  e.g.,  a  b  =  a1  gift . 


2.3.  Thermodynamic  framework 

We  assume  a  Helmholtz  free  energy  potential  of  the  form 

^  =  ^(ee,0,a  (10) 

where  the  intermediate  configuration  elastic  strain  2(ee)ajS  =  fTSabff  ~  with  g ^  a  metric  tensor  on  b , 
which  in  practice  is  chosen  as  Kronecker’s  delta  for  simplicity  (Simo  and  Ortiz,  1985),  a  typical  assump¬ 
tion  in  finite  elastoplasticity  theory  (Clayton  et  al.,  2004).  The  symbol  £  denotes  a  dimensionless  scalar 
internal  variable  representing  stored  micro-elastic  energy  associated  with  crystal  defects  that  may  impede 
shearing  on  each  slip  system.  Expanding  the  time  rate  of  f  and  inserting  the  result  into  inequality  (9)  yields 

tr(<7g(r  + 19  +  P))  -  >  p(0£^ee  +  (0^)0  +  (0^K  +  ft/),  (11) 

where  the  scalar  product  of  two  second  rank  objects  is  defined  by  A:B  =  A^By,  and  where  the  subscript  fol¬ 
lowing  the  6-operator  denotes  partial  differentiation  with  respect  to  the  subscripted  variable.  Additional 
algebraic  manipulations  give 

?  (se  -  p8<-i/0  :  ee  +  (tr(ffg)aT  -  p(de  «A)  -  pn)& ~  p(d^)t  +  ^  ,  (12) 

P  a=l  U 

where  p  is  the  mass  density  in  configuration  b  and  the  elastic  second  Piola-Kirchhoff  stress 
(se)^  =  ff^-laaabfeb~ip  with  f  =  p/p  =  detfe.  The  resolved  Cauchy  stress  on  slip  system  a  is  found  by 
T(a)  =  <r:(gs(a)(g)m(a)).  Assuming  that  the  elastic  strain  ee  and  temperature  9  can  be  varied  apart  from  the 
inelastic  rates  y^  and  \  (cf.  Coleman  and  Noll,  1963),  (12)  is  satisfied  when 


Se  =  pdee\J/, 

(13) 

Qfy 

n  =  —  tr(ffg)  -del y 

(14) 

=x 

a=l 

(15) 

Notice  that  the  first  term  on  the  right  of  (14),  denoted  by  % ,  arises  as  a  consequence  of  the  explicit  inclusion 
of  thermal  deformation  in  the  kinematic  description  (1).  Rearranging  the  energy  balance  in  (6)  by  appealing 
to  (13)  and  (14)  gives 

n 

pOf  —  pOdei/i  =  pr  —  divq  +  ^  f^y^  —  p(0^)£. 

a=l 


(16) 
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The  specific  heat  capacity  c  is  defined  as 

c  =  dee=  dtle(dgri)  =  -9,  (dee^)  +  ^edgX  (17) 

where  we  have  invoked  the  differential  of  relation  (8).  Assuming  isotropic  heat  conduction  in  the  spatial 
frame  dictated  by  Fourier’s  law  (Zhou  et  ah,  1994) 

q  =  -kVx6  (18) 

with  k(x,0)  the  spatial  thermal  conductivity,  we  can  rewrite  the  localized  energy  balance  (16)  as 

pcO  =  f>(a)y(a)  -p((0^)  -  0(doz'l'))t+  p0d6eeil/:ee  +  di v(kVx0)  +  pr  .  (19) 

si—  s  f  ^  v  _  y  s  y  v  y  s-_l  s 

^  a=l  ^  - - -  v'  ^ 

temperature  change  s ^  energy  of  lattice  defects  thermoelastic  coupling  heat  conduction  heat  supply 

plastic  dissipation 

2.4.  Constitutive  theory 

A  free  energy  potential  per  unit  intermediate  configuration  volume  is  specified  as 

W  =  \k0(9)92  -  l-K,y  +  n(P)c?  :  ee'  +y(9),  (20) 

where  =  tr(ee)  measures  the  elastic  volume  change,  K0,  Ki,  and  p  are  elastic  stiffness  constants  and  k  is  a 
dimensionless,  material-dependent  scalar  parameter  that  we  assume  is  independent  of  strain  rate  and  tem¬ 
perature.  We  shall  later  demonstrate  the  explicit  relationship  between  k  and  the  fraction  of  plastic  dissipa¬ 
tion  converted  to  heat  energy  (Eq.  (37)).  The  covariant  deviatoric  elastic  strain  tensor  in  the  intermediate 
frame  is  defined  by  ee'  =  ee  —  (l/3)$g.  The  function  y(0)  =  —c6ln(6/60)  accounts  for  the  purely  thermal 
energy  (Rosakis  et  al.,  2000),  with  0o  a  reference  temperature  at  which  y  =  0.  From  Eq.  (13)  and  partial 
differentiation  of  (20),  we  see  that  the  stress  tensor  satisfies 

se  =  =  (Ko-Krf)  9g~l  +2pg~lee'g~\  (21) 

=K 

where  K  is  the  temperature-  and  volume-dependent  (i.e.,  pressure-dependent)  effective  bulk  modulus.  The 
hydrostatic  pressure  p  and  deviatoric  stress  se  in  the  intermediate  frame  become,  from  (21), 

p  =  - 1  tr(se)  =  -K9,  se'  =  se  +pg~'  =  2/ig_1eeg_1.  (22) 

Thermoelastic  properties  for  W  and  binder  phases  are  listed  in  Table  1,  where  the  temperature  change  is 
defined  by  A6  =  6  —  273  K.  The  dependencies  of  the  elastic  moduli,  conductivity,  specific  heat,  and  coeffi- 


Table  1 

Thermoelastic  properties 


Property 

Value  (W) 

Value  (matrix) 

li 

160.3  -  O.O146(A0)  -  O.OOOOO33(A0)2  GPa 

99.0  GPa 

K0 

310.8  -  O.O1759(A0)  -  O.OOOOOO3(A0)2  GPa 

190.3  GPa 

Ki 

1450  GPa 

- 

GC’j' 

5.3(10)  6  +  1.584(10r9  (A0)/K 

1.5(10r5/K 

C 

135.8(1  -  48O5/02)  +  0.009120  +  0.0023 1(  1O)_603  J/(kg  K) 

382  J/(kg  K) 

k 

160  -  O.O5(A0)  W/(m  K) 

1 00  W/(m  K) 

Po 

19350  kg/m3 

9200  kg/m3 
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cient  of  thermal  expansion  of  the  pure  W  phase  on  temperature  were  obtained  from  Yih  and  Wang  (1979), 
most  valid  for  0  ^  1500  K.  Temperature  and  pressure  dependencies  of  thermoelastic  properties  of  the  ma¬ 
trix  phase  are  unknown  and  were  neglected  in  the  computations.  The  value  of  the  bulk  modulus  term  K\ 
was  obtained  from  a  linear  fit  to  the  data  of  Ruoff  et  al.  (1998),  most  valid  for  applied  pressures  less  than  35 
GPa.  In  the  final  row  of  Table  1,  the  mass  density  in  the  reference  configuration  is  written  as  po,  with  the 
listed  value  corresponding  to  a  reference  temperature  of  6  =  300  K. 

Note  that  single  crystalline  W  is  unusual  in  the  sense  that  it  is  elastically  isotropic  at  the  reference  tem¬ 
perature  and  low  pressures.  Presumably,  single  crystals  of  the  FCC  binder  phase  are  not  elastically  isotro¬ 
pic;  however,  we  resort  here  to  isotropic  stiffness  constants  obtained  from  elastic  wave  speed  measurements 
on  poly  crystalline  sample  (Zhou,  1993),  as  anisotropic  stiffness  constants  are  presently  unavailable  for  sin¬ 
gle  crystals  of  the  binder  phase.  Note  also  that  pressure  dependency  of  the  shear  moduli  is  neglected  here,  as 
this  effect  is  relatively  unimportant  in  the  present  context  where  plasticity  dominates  the  deviatoric  stress- 
strain  response.  On  the  other  hand,  a  volume-dependent  bulk  modulus  is  necessary  in  order  to  numerically 
capture  the  sharpening  of  elastic  stress  waves  into  a  shock  discontinuity  commensurate  with  the  increase  in 
wave  speed  occurring  in  highly  compressed  material  behind  a  pressure  front  (Becker,  2004).  Our  implemen¬ 
tation  is  relatively  simple,  involving  only  a  linear  dependence  of  bulk  modulus  on  elastic  volume  change  and 
a  quadratic  dependency  of  stiffness  upon  temperature.  This  approach  is  deemed  sufficient  over  the  present 
range  of  simulation  conditions  (p  <  0.1  K0).  More  robust  formulations  involving  equations-of-state  (cf. 
Griineisen,  1926)  and  pressure  derivatives  of  shear  moduli  may  be  necessary  for  modeling  elastically  aniso¬ 
tropic  crystals  deformed  at  pressures  on  the  order  of  the  reference  bulk  modulus — see  for  example  simula¬ 
tions  of  poly  crystalline  tantalum  by  Becker  (2004).  Note  also  that  the  stiffness  constants  and  elastic  volume 
change  ft  are  defined  with  respect  to  the  intermediate  configuration,  a  conventional  practice  in  finite  crystal 
plasticity  theory  (Asaro,  1983).  This  poses  no  conceptual  difficulties  and  is  rigorous  if  the  material  param¬ 
eters  are  determined  in  a  manner  consistent  with  their  definitions.  However,  in  experiments  the  constants 
are  often  determined  by  perturbing  the  material  from  its  equilibrium  state  in  the  current  configuration.  Spe¬ 
cifically,  data  used  to  determine  the  volume-dependent  parameter  K\  (Ruoff  et  al.,  1998)  were  given  in  terms 
of  Cauchy  pressure  p  =  — (1/3)  tr  and  elastic  volume  change  per  unit  reference  volume  in  the  absence  of 
plasticity  or  thermal  expansion,  /  -  1.  The  final  value  of  Kx  appearing  in  Table  1  was  found  from  a  con¬ 
version  using  the  relations  ft  =  1.017(/  -  1)  and  p  =  1.036 p  when  f  =  0.9,  corresponding  to  a  10%  volu¬ 
metric  reduction  at  an  applied  Cauchy  pressure  of  36.6  GPa.  The  linear  elasticity  constants  (K0  and  p) 
of  Table  1  correspond  to  those  measured  at  small  elastic  strains,  i.e.  when  /  «  1. 

Our  decomposition  of  the  deformation  gradient  into  distinct  elastic  and  thermal  terms  differs  somewhat 
from  usual  finite  thermoelasticity  (cf.  Marsden  and  Hughes,  1983).  Consider  a  purely  thermoelastic  process 
at  a  material  point  where  the  total  deformation  is  prohibited,  such  that  f  =  f  f°  =  1,  and  the  temperature  is 
increased  above  its  reference  value.  Then  due  to  thermal  expansion,  detf0  >  1,  and  in  order  to  compensate 
and  maintain  null  total  strain,  detF  <  1  and  hence  ft  <  0,  resulting  in  a  positive  pressure  p  >  0  from  (22). 

Plasticity  models  are  discussed  next.  A  power-law  viscoplastic  flow  rule  (Hutchinson,  1976)  is  chosen  to 
model  the  time  rate  of  plastic  deformation  within  each  phase  of  the  bulk  material: 

*  sgn(TW).  (23) 

Here  y0  and  m  are  material  parameters,  g (a)  is  the  slip  resistance  due  to  dislocation  barriers,  =  f  is 
the  projected  shear  stress  pulled  back  to  the  intermediate  configuration  b  of  Fig.  1,  and  sgn(x)  =  x/\x\,  with 
sgn(0)  =  1.  Thermal  softening  attributed  to  increased  dislocation  mobility  at  high  temperatures  is  incorpo¬ 
rated  via  the  power-law  form  (Klopp  et  al.,  1985) 


=g{o\e/e0y 


(24) 
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with  g Qa)  the  flow  resistance  at  reference  temperature  60  (a  material  parameter)  and  p  a  dimensionless  con¬ 
stant  for  each  material.  We  postulate  the  following  relationship  between  the  “average”  hardening  over  all 
systems  at  fixed  reference  temperature  and  the  internal  variable 

l-T,(g{l)-g)))  =  ^(00)b^,  (25) 

a=1 

with  g ^  an  initial  yield  stress,  b  the  magnitude  of  the  Burgers  vector  in  the  reference  lattice,  and  pT  the  total 
dislocation  line  length  per  unit  intermediate  configuration  volume  associated  with  shearing  impedance.  The 
square-root  dependence  of  flow  stress  on  dislocation  density  is  a  typical  assumption  in  the  plasticity  and 
materials  science  literature  (Taylor,  1934;  Kuhlmann-Wilsdorf,  1985;  Zikry  and  Kao,  1996),  as  is  the  asser¬ 
tion  of  linear  dependence  of  stored  lattice  energy  on  dislocation  density  (Bammann,  2001;  Svendsen,  2002) 
implied  jointly  by  Eqs.  (20)  and  (25).  The  scalar  proportionality  factor  a  accounts  for  dislocation  interac¬ 
tions  (Kobytev  et  al.,  1984;  Ashmawi  and  Zikry,  2003).  Both  lattice  friction  stress  (Qiu  et  al.,  2001)  and 
effects  of  initial  dislocation  density  are  incorporated  in  the  initial  yield  stress  g^\ 

Specifically  for  pure  W  crystals,  we  allow  slip  in  the  (111)  close-packed  directions  on  any  of  the  {110} 
and  {112}  families  of  planes,  meaning  the  number  of  potentially  active  slip  systems  is  n  =  24  (Subhash 
et  al.,  1994).  Possible  slip  on  {123}  families  of  planes,  thought  to  be  inactive  at  room  temperature  (Argon 
and  Maloof,  1966;  Subhash  et  al.,  1994),  is  not  represented  in  our  model.  Evolution  of  slip  resistance  at 
reference  temperature  60  is  dictated  by  a  hardening-minus-dynamic-recovery  relation  reminiscent  of  a  form 
introduced  in  a  macroscopic  plasticity  setting  by  Armstrong  and  Frederick  (1966)  and  more  recently  used  in 
a  crystal  plasticity  framework  by  Horstemeyer  et  al.  (1999): 

go]  =  A  qxp  |  yw  |  -Bgf  I  V®  |,  (26) 

n=i  H=1 

with  the  interaction  matrix  satisfying 

^*^  +  *(1-^),  (27) 

where  q  is  the  latent  hardening  ratio.  Material  parameters  for  the  plasticity  models  are  compiled  in  Table  2. 
As  calibrations  to  experimental  data  used  to  obtain  these  parameters  are  discussed  in  detail  in  a  previous 
publication  by  the  author  (Clayton,  2005),  only  a  brief  overview  is  given  here.  Values  of  g^  =  gg  |,=0,  A, 
and  B  in  Eq.  (26)  were  chosen  to  fit  the  quasi-static  compression  data  for  polycrystalline  W  of  Dummer 
et  al.  (1998),  and  a  relatively  large  value  of  q  =  1.4  was  selected  for  implementation  in  Eq.  (27)  following 
the  discussion  on  strong  latent  hardening  in  W  single  crystals  by  Horwath  (1994).  Strain  rate  and  thermal 


Table  2 

Properties  for  inelastic  material  response 


Parameter 

Value  (W) 

Value  (matrix) 

7o 

0.001 

0.001 

m 

20 

20 

q 

1.4 

1.0 

A 

630  MPa 

200  MPa 

B 

1.5 

0.4 

4a) 

500  MPa 

150  MPa 

p 

-1.5 

-1.5 

®0 

300  K 

300  K 

a 

0.439 

1.03 

b 

0.275  nm 

0.257  nm 

K 

1333 

200 
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sensitivity  parameters,  m  and  p  respectively,  were  obtained  from  the  crystal  plasticity  model  for  pure  W  of 
Lee  et  al.  (1999).  The  value  of  the  reference  strain  rate  y0,  required  on  dimensional  grounds  (Hutchinson, 
1976),  was  chosen  for  convenience  during  the  process  of  selecting  the  other  parameters.  Aggregates  of  300 
initially  randomly-oriented  crystals  subjected  to  uniform  deformation  constraints  (Taylor,  1938)  in  the  ab¬ 
sence  of  inertia  were  used  in  the  calibrations. 

In  the  W-Ni-Fe  matrix  material,  the  number  of  potentially  active  slip  systems  was  chosen  as  n  =  12.  We 
assume  that  dislocations  glide  in  (1 1 0)  close-packed  directions  on  { 1 1 1 }  planes  for  this  FCC  metal.  Strain 
rate-  and  temperature-dependent  slip  resistances  in  crystals  comprising  the  more  compliant  binder  phase 
evolve  via  relations  (26)  and  (27),  though  with  different  values  of  A,  B,  and  q  than  those  invoked  for  the 
pure  W  grains.  Specifically,  constants  A  and  B  in  (26)  and  initial  conditions  g ^  were  chosen  to 

fit  the  quasi-static  (y  =  10-4)  torsion  data  of  Zhou  (1993),  and  the  strain  rate  sensitivity  exponent  was  as¬ 
sumed  to  be  the  same  as  that  for  the  pure  W  grains  following  the  experimental  data  of  Zhou  (1993).  We  also 
selected  y0  to  match  the  value  for  the  pure  W,  and  used  Taylor’s  (1934)  latent  hardening  assumption  typical 
for  FCC  metals,  q  =  1.0.  Choices  of  a  and  k  for  each  phase  are  discussed  below  in  Section  2.5. 

2.5.  Implementation  and  internal  variables 


An  implicit  hyperelastic-viscoplastic  algorithm  (cf.  Cuitino  and  Ortiz,  1992)  was  employed  to  integrate 
the  elastic-plastic  constitutive  response.  Let  subscripts  t  and  t  +  At  denote  consecutive  computation  times 
in  a  nonlinear  analysis,  i.e.,  start  and  end  times  in  a  particular  iteration  cycle.  Slip  rates  for  a  given  incre¬ 
ment  spanning  times  t  and  t  +  At  are  found  implicitly  from  values  of  the  resolved  shear  stress  and  harden¬ 
ing  variables  at  the  end  of  the  cycle: 


y(a)  =  To 


=(«) 

lt+At 


(a) 

St+At 


(28) 


We  use  an  iterative  procedure  to  solve  (28),  since  and  g^At  depend  upon  the  solution  variables  y^\ 
Notice  that  T^At  and  gf^At  depend  upon  6,  via  Eq.  (24)  and  the  temperature  dependence  of  elastic  moduli. 
The  temperature  rate  for  a  given  increment  spanning  t  and  t  +  Ans  found,  from  Eq.  (19)  with  r  =  0,  explic¬ 
itly  in  terms  of  quantities  at  time  t : 


6  =  1 4  V  T(<,)y(a)  +  ?0„(se  :  ee)  +  L div(JfcVx0) 
\pc  jrt  c 


pc 


where 


P=\Y  T(x)y(x)  -  p((drf)  -  0(de(m)  ( f)  ^ 
V  a=l  I  V  a=l 


-1 


(29) 


(30) 


meaning  that  1  -  /?  is  the  ratio  of  the  rate  of  stored  energy  of  lattice  defects  to  the  rate  of  plastic  work 
(Taylor  and  Quinney,  1934).  In  practice,  /?  is  often  prescribed  a  fixed  magnitude  for  a  particular  material 
regardless  of  the  temperature-deformation  history,  with  0.8  ^  /?  ^  1.0  typical  for  engineering  metals  (Ara- 
vas  et  al.,  1990;  Zhou  et  al.,  1994).  We  do  not  employ  this  constitutive  assumption  here.  The  temperature  at 
time  t  +  At  is  updated  simply  via 


0t+At  —  0t  +  OAt , 


(31) 


and  the  thermal  deformation  gradient  at  the  end  of  the  step  is  found  from 

Cm  =  exp(aT01Af)f?, 


(32) 
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with  exp()  the  matrix  exponential  function  and  the  property  aT  in  (32)  evaluated  at  0t+At •  For  a  given  time 
increment,  Eqs.  (28)  are  solved  implicitly  using  values  of  6 ,  9t+ At,  and  fet+At  found  from  (29)-(32).  The  ther¬ 
moelastic  term  in  (19),  upon  assuming  f  =  detf0  «1  +  3aT(0  —  0O)  can  be  rewritten,  to  first-order,  as 

pOdo#il/  :  ee  =  0f~%( se  :  ee)  -  p6j9de(je-l)(se  :  ee) 

»  Of-'  { [Mr^  +  M(rxgps + i'V*)]  (*%  +  i33) 

where  for  pure  W, 

=  -3.4  +  O.OO65(A0)  MPa/K,  00iu  =  -14.6  -  O.OO66(A0)  MPa/K,  (34) 

and  where  the  second  of  (34)  agrees  with  the  form  given  in  Table  1.  Temperature  dependencies  of  elastic 
moduli  of  the  matrix  phase  are  unknown  and  hence  were  neglected  in  the  computation  of  thermoelastic 
heating  effects. 

Additional  explanation  regarding  the  role  of  the  internal  variable  state  variable  in  each  phase  corre¬ 
sponding  to  the  total  dislocation  density,  is  now  merited.  Manipulating  Eq.  (25)  and  then  differentiating 
with  respect  to  time  give,  respectively, 


S3 vs  =  — =  ^- i>?’ - 4”).  5  = 

Kfi  a pn  “  7  ^\[Py 

Additionally, 

POM)  =  Odeiwb^)  =  K6(den)(by/pd  =  g  (g^-gf)- 

Substituting  (35)  and  (36)  into  (30)  then  yields 


0=1- 


ify  +  8(defj)) 

(oifiit)2 


a=l 


#=1 


(«) 


a=l 


(35) 


(36) 


(37) 


The  dislocation  density  variable  can  be  determined  from  (35)  once  a  is  known,  i.e., 


Pt  — 


5c[inb 


(38) 


Notice  that  (38)  is  a  convenient  relationship  between  the  total  dislocation  density  pT  and  the  hardness 
averaged  over  n  potential  slip  systems  at  reference  temperature  60.  The  value  of  parameter  a  was  chosen 
upon  consideration  of  published  data  on  the  dislocation  line  density  (Argon  and  Maloof,  1966).  It  should 
be  noted  that  since  the  experimental  data  we  used  to  justify  our  choice  of  a  suffers  from  an  acknowledged 
lack  of  precision  (Argon  and  Maloof,  1966),  our  variable  pT  should  be  viewed  as  a  qualitative,  yet  physi¬ 
cally-based,  measure  useful  for  comparing  the  relative  degree  of  strain  hardening  accumulated  in  regions  of 
W  and  matrix  phases.  Furthermore,  since  the  parameter  a  was  chosen  from  calibration  of  pT  with  exper¬ 
imental  dislocation  density  measurements  (see  Clayton,  2005),  the  variable  pT  implicitly  enters  the  energy 
balance  through  the  presence  of  a  in  Eq.  (37). 

Caution  must  be  taken  when  selecting  a  value  of  k  in  order  to  ensure  satisfaction  of  the  reduced  dissi¬ 
pation  inequality  (15).  If  data  revealing  evolution  of  /?  is  available  from  physical  experiments  (Hodowany 
et  al.,  2000;  Rosakis  et  al.,  2000),  then  Eq.  (37)  can  be  inverted  to  select  an  appropriate  value  of  k. 
Vice-versa,  if  k  is  known,  then  instantaneous  values  of  /?  can  be  predicted  via  (37).  From  (20)  and  (37) 
we  see  that  /?  — ►  1  as  k  — ►  0,  meaning  no  energy  of  cold  working  is  stored  in  the  lattice  when  k  —  0.  Subhash 
et  al.  (1994)  deduced  from  infrared  temperature  measurements  that  /?  approached  unity  in  high-rate 
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(i.e.,  Kolsky  bar)  compression  tests  on  pure  W  polycrystals,  while  Zhou  et  al.  (1994);  Zhou  (1998)  invoked 
P  =  0.90  for  binder  and  pure  W  phases  in  dynamic  finite  element  calculations.  Notice  that  values  of  /? 
depend  strongly  upon  the  form  of  the  rate  equations  for  the  slip  resistances  g^\  the  latter  which,  upon 
consideration  of  (38),  can  be  interpreted  as  evolution  laws  for  the  dislocation  density.  Additional  detail 
could  be  injected  into  (20)  to  reflect  complex  mechanisms  of  energy  storage,  but  more  experimental  data 
are  needed  to  justify  a  more  advanced  treatment. 


3.  Fracture  modeling 


Consistent  with  experimental  observations  of  fracture  in  tungsten  alloys  (Zamora  et  al.,  1992;  Wee- 
rasooriya  et  al.,  1994;  O’Donnell  and  Woodward,  2000;  Woodward  and  O’Donnell,  2000;  Weerasooriya, 
2003),  in  our  framework  new  surfaces  may  be  created  at  interfaces  between  grains  of  like  or  different 
phases,  i.e.  along  W-W  grain  boundaries  and  W-matrix  boundaries.  Intragranular  cleavage  and  separation 
at  matrix-matrix  grain  boundaries  are  not  represented,  as  these  phenomena  were  observed  less  frequently 
in  the  aforementioned  high-rate  tests  (Weerasooriya  et  al.,  1994;  Weerasooriya,  2003)  and  in  the  room  tem¬ 
perature,  quasi-static  experiments  of  Woodward  and  O’Donnell  (2000).  We  invoke  the  dynamic  finite  ele¬ 
ment  approach  (see  e.g.  Ortiz  and  Pandolfi,  1999,  for  theoretical  underpinnings),  with  new  fracture  surfaces 
emerging  at  interfaces  between  continuum  elements  when  the  traction  resolved  on  the  potential  initiation 
site  exceeds  the  intrinsic  strength  of  the  interface  (here,  a  temperature-dependent  material  parameter). 
Hence,  fracture  initiates  when  one  of  the  following  local  stress-based  criteria  is  attained: 

s  =  So  +Si(0  -  0o),  t  =  T0  +  ti(0  -  0O),  (39) 


where  £  and  t  are  the  resolved  normal  traction  and  shear  traction  on  the  interfacial  surface,  measured  per 
unit  reference  configuration  area,  and  sq  and  fo  are  material  parameters  specifying  the  normal  and  tangen¬ 
tial  (i.e.,  mode  I  and  mode  II)  fracture  strengths  of  the  interface  at  reference  temperature  60.  Linearly  tem¬ 
perature-dependent  initiation  strengths — capturing  ductile-brittle  transitions,  for  example — are  specified 
by  prescribing  nonzero  values  of  Si  and  (dimensions  of  stress/temperature),  and  various  mixed-mode  ini¬ 
tiation  criteria  are  readily  enforced  by  changing  the  ratio  S0/t0.  In  the  finite  element  implementation,  dupli¬ 
cate  nodes  are  generated  along  all  potential  fracture  surfaces  during  the  meshing  stage.  Initially  coincident 
nodes  share  the  same  velocity  and  temperature  histories  until  either  of  conditions  (39)  is  reached. 

The  response  of  the  degraded  material  at  interfaces  following  initiation  is  specified  by  the  following  cou¬ 
pled  irreversible  traction-displacement  laws  (see  e.g.  Camacho  and  Ortiz,  1996): 


s  =  (so+.si(0-0o))^l 
s  =  (s 0  +  S i(0  —  0o))  ^1 

T  =  (tO+Ti(0-0O))^1 
T=  (tO+Ti(0—  0O))(l 


j-'j  >  (loading) 

T1)  7T  ’  (unloading) 
J  Onl 


* 

|*i 

* 


1  _^)sgn(<5t)’ 

(&n)\  * 

v  *Mii’ 


(loading) 

(unloading) 


(40) 


(41) 


with  3n  and  3t  the  normal  and  tangential  crack  opening  displacements,  <5nl  and  3n  the  maximal  values  of  3n 
and  3t  achieved  during  the  loading  history,  and  3C  a  material  parameter  describing  the  separation  distance 
beyond  which  the  cohesive  interface  no  longer  supports  traction  (i.e.  complete  separation).  The  Macaulay 
bracket  is  written  (),  satisfying  (x)=xVx^0  and  (x)  =  0  Vx  <  0.  As  shown  in  Fig.  2(a)  for  pure  tension 
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(a)  (b) 

Fig.  2.  Isothermal  cohesive  laws  for  pure  tension  (a)  and  pure  shear  (b). 


and  Fig.  2(b)  for  pure  shear,  unloading  (reloading)  occurs  linearly  to  (from)  the  origin  in  traction-displace¬ 
ment  space.  Also  prior  to  attainment  of  either  of  criteria  (39),  heat  conduction  proceeds  as  usual.  Upon 
initiation  of  damage  and  subsequent  separation,  however,  we  enforce  the  null  flux  conditions  Vx0- n  =  0, 
with  n  the  outward  normal  to  the  newly-created  fracture  surface. 


4.  Simulations  and  results 

The  constitutive  models  presented  in  Sections  2  and  3  for  bulk  and  interfacial  responses,  respectively, 
were  implemented  within  the  EPIC  Lagrangian  finite  element  solver  (Johnson  et  al.,  1997,  2001;  Schoenfeld, 
1998).  In  this  approach,  the  equations  of  motion  are  integrated  directly  and  explicitly  (see  for  example  Bely- 
tschko  et  al.  (1976)).  The  deformation  gradient  is  updated  within  each  integration  point  as 

ft+At  =  exp(lA?)f,,  (42) 

where  1  =  ff-1  is  the  velocity  gradient  that  is  assumed  constant  over  the  time  interval  (t,  t  +  At).  The  con- 
stitutive  update  proceeds  at  each  integration  point  via  the  methodology  outlined  in  Section  2.5.  Contribu¬ 
tions  to  the  energy  balance  (19)  from  plastic  dissipation,  lattice  defects,  and  thermoelastic  coupling  are 
calculated,  and  then  the  temperature  field  is  updated  explicitly,  following  the  procedure  of  Johnson 
(1981).  Possible  contributions  from  cohesive  elements  to  global  mechanical  and  thermal  force  vectors  are 
then  calculated  prior  to  the  enforcement  of  external  boundary  conditions  and  initiation  of  the  integration 
step  for  the  deformation  field  of  the  next  cycle.  The  stable  time  increment  for  each  cycle  is  chosen  as  a  small 
fraction  of  the  time  needed  for  a  longitudinal  elastic  stress  wave  to  traverse  the  smallest  element  in  the 
mesh.  A  more  detailed  description  of  the  explicit  finite  element  implementation  is  given  in  Appendix  A. 

Two-dimensional  calculations  were  conducted,  with  micro  structures  for  meshing  obtained  from  cross- 
sectional  samples  of  a  kinetic  energy  (KE)  rod  of  the  WHA  material.  Fields-of-view  from  two  different 
microstructures  were  considered,  as  shown  in  Figs.  3  and  4,  each  a  square  window  of  size  L  =  150  pm. 
As  a  result  of  processing,  W  grains  in  the  longitudinal  section  are  elongated  in  a  preferred  direction  relative 
to  those  in  the  transverse  section,  with  grains  in  the  latter  appearing  nearly  equiaxed.  Lattice  orientations 
used  as  initial  conditions  in  some  simulations  were  obtained  from  Electron  BackScatter  Diffraction  (EBSD) 
techniques  (Humphreys,  2001).  In  Figs.  3(b)  and  4(b),  various  colors  indicate  orientations  of  normal  vec¬ 
tors  to  the  crystallographic  planes  exposed  at  the  surface  (for  clarity,  only  the  W  phase  is  shown  here).  The 
meshes  consist  of  constant  strain  triangular  elements,  generated  with  the  PPM200F  software  (Langer 
et  al.,  2003),  featuring  algorithms  for  selective  grid  refinement  along  material  interfaces.  The  mesh  of 
Fig.  3(c)  features  52  distinct  W  grains,  a  volume  fraction  of  82.1%  pure  W,  16999  standard  3-node  bulk 
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Fig.  3.  Longitudinal  section:  (a)  penetrator,  (b)  lattice  orientation  map,  (c)  FE  mesh. 


Fig.  4.  Transverse  section:  (a)  penetrator,  (b)  lattice  orientation  map,  (c)  FE  mesh. 


elements,  and  up  to  1879  4-node  cohesive  elements  inserted  along  grain  and  phase  boundaries.  The  grid 
shown  in  Fig.  4(c)  features  51  W  grains,  83.9%  pure  W,  16834  bulk  elements,  and  up  to  2028  cohesive 
elements. 

The  simulations  reported  here  describe  the  behavior  of  the  material  under  impulsive  compression  load¬ 
ing,  similar  to  conditions  that  would  occur  in  high  velocity  plate  impact  tests  (Dandekar  and  Weisgerber, 
1999).  Our  intention  is  to  track  the  effects  of  the  bulk  micro  structural  and  interfacial  properties  on  the  local 
thermomechanical  response  history,  effects  which  are  not  easily  deduced  from  macroscopic  physical  exper¬ 
iments  alone.  Via  appropriate  choice  of  boundary  conditions,  we  are  able  to  create  a  situation  in  which  the 
material  spalls  as  a  result  of  wave  reflection,  a  setting  in  which  we  can  compare  the  fracture  behavior  with 
other  spall  experiments  of  the  WHA  material  reported  in  the  literature  (Baoping  et  al.,  1994;  Dandekar  and 
Weisgerber,  1999)  and  parametrically  investigate  the  effects  of  local  interfacial  fracture  properties  on  the 
response  (Espinosa  and  Zavattieri,  2003a, b). 

The  following  boundary  and  initial  conditions  were  specified.  Analyses  were  plane  strain  in  the  xl  -  x2 
plane,  in  essence  representing  columnar  polycrystals  extended  infinitely  in  the  out-of-plane  direction.  How¬ 
ever,  out-of-plane  elastic  and  plastic  deformations  were  permitted  (i.e.,  the  3D  material  models  discussed  in 
Section  2  were  employed,  with  the  full  complement  of  slip  systems  enabling  out-of-plane  lattice  rotations), 
so  long  as  the  total  deformation  field  remained  planar.  We  denote  the  lower  and  upper  edges  of  the  domain 
by  X2  =  0  and  X2  =  L,  respectively.  And  let  the  left  and  right  edges  be  denoted  by  X1  =  0  and  X1  =  L, 
respectively.  The  velocity  initial  and  boundary  conditions  are  listed  as 
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X 


1 


(250  m/s)(t/l(T9  s) 


for  0  ^  t  ^  10  9  s 
along  X1  =  0 


250  m/s<^ 


for  10  9s  ^  ^  10  8  s 
along  X1  =  0 


(250  m/s)([l.l(10  8)  s  —  ^] / 10  9  s) 

0  f  for  t  ^  1 . 1  ( 10-8)  s 
\  along  X1  =  0. 


for  10“8  1.1(10“8)  s 

along  X1  =  0 


(43) 


Conditions  (43)  result  in  a  rightward-moving  compressive  pulse  of  duration  9  ns.  Additionally,  boundaries 
along  X2  =  0  and  X2  =  L  were  free  of  shear  stress  but  fixed  in  the  x2-direction,  thus  establishing  uniaxial 
strain  conditions  similar  to  those  that  would  be  encountered  by  a  small  volume  element  of  the  crystalline 
material  embedded  within  the  target  specimen  of  a  plate  impact  test  (see  e.g.  Meyers,  1994;  Becker,  2004): 


tl  =  onri2  =  0 
.  2  =  0  along 


(x2  =  o 

\x2  =  l. 


(44) 


The  right-hand  side  of  the  mesh  was  traction  free,  i.e.  a'J  tij  =  0  along  Xx  =  L  for  t  ^  0.  Null  heat  flux 
boundary  conditions  were  enforced  as 

(X1  =0 


Vx6  •  n  =  0  along  < 


Xl=L 
X2  =  0 

lx2=z, 


(45) 


with  n  the  outward  normal  vector  in  the  spatial  frame,  and  an  initial  temperature  60  =  300  K  was  specified 
uniformly  throughout  the  domain. 

Results  from  nine  distinct  compression  simulations,  denoted  as  simulations  1-9  in  Table  3,  are  reported 
here.  Micro  structures  and  initial  lattice  orientations  are  labeled  in  columns  2  and  3  of  Table  3.  Simulations 
1-3  and  6-9  invoke  as  initial  conditions  the  experimentally-obtained  crystallographic  textures  correspond¬ 
ing  to  Figs.  3(b)  and  4(b),  while  simulations  4  and  5  feature  purely  random  initial  textures  for  the  W  grains 
and  uniform  orientations  for  the  matrix  phase.  Interfacial  fracture  properties  were  also  varied  among  sim¬ 
ulations,  as  indicated  in  columns  4-6.  For  simulations  in  which  fractures  were  permitted,  values  of 


Table  3 
Simulations 


Simulation 

Microstructure 

Lattice  orientation 

£o  (GPa) 

to  (GPa) 

Si  =  ii  (GPa/K) 

1 

Longitudinal 

Experimental  1 

2.0 

oo 

0 

2 

Longitudinal 

Experimental  1 

2.0 

4.0 

0 

3 

Transverse 

Experimental  2 

2.0 

oo 

0 

4 

Longitudinal 

Random  1 

2.0 

oo 

0 

5 

Transverse 

Random  2 

2.0 

oo 

0 

6 

Longitudinal 

Experimental  1 

2.0 

oo 

0.01 

7 

Transverse 

Experimental  2 

2.0 

oo 

0.01 

8 

Longitudinal 

Experimental  1 

oo 

oo 

0 

9 

Transverse 

Experimental  2 

2.0 

2.0 

0 
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so  =  2.0  GPa  were  nominally  assigned,  in  order  to  roughly  match  peak  stress  levels  reported  in  experimental 
high-rate  macroscopic  data  of  Weerasooriya  (2003).  Dandekar  and  Weisgerber  (1999)  estimated  from  plate 
impact  experiments  the  spall  threshold  stress  of  WHA  to  he  between  1.7  and  2.0  GPa.  Stochastic  variations 
in  strength  and  toughness  among  interfaces  are  expected  (Woodward  and  O’Donnell,  2000);  though  these 
aspects  are  not  addressed  in  the  present  work,  such  variations  could  be  readily  implemented  within  a  cohe¬ 
sive  finite  element  framework  such  as  ours  (cf.  Zhou  and  Molinari,  2004).  As  indicated  in  column  5,  mixed¬ 
mode  initiation  criteria  were  investigated  in  several  cases  by  prescribing  various  values  of  the  initiation 
strength  t0.  For  simulations  in  which  shear  stress-induced  initiation  was  prohibited  (i.e.,  f0  of  oo),  the  value 
of  To  in  (41)  was  set  equal  to  So  immediately  after  the  normal  stress  initiation  criterion  was  achieved,  in  order 
to  avoid  infinite  internal  reaction  forces  due  to  shearing  resistance.  The  choice  of  infinite  initial  shear 
strength  is  an  idealization  motivated  by  Kolsky  bar  experiments  (Weerasooriya  and  Beaulieu,  1993;  Wee¬ 
rasooriya  and  Moy,  1998;  Weerasooriya,  2003)  that  indicate  a  tendency  for  mode  I  fracture  to  occur  readily 
relative  to  modes  II  and  III  in  the  WHA  material  system.  Thermally-dependent  fracture  strengths  were 
investigated  in  simulations  6  and  7,  in  which  we  allowed  the  fracture  strengths  to  increase  linearly  with 
increasing  temperature,  in  association  with  the  brittle-to-ductile  transition  in  tungsten.  Such  increase  is 
motivated  by  the  experimental  observation  of  reduced  frequency  of  intergranular  fracture  events  at  higher 
ambient  temperatures  (O’Donnell  and  Woodward,  2000)  with  our  designation  of  a  linear  dependency  of 
strength  on  temperature,  and  particular  value  of  Si,  following  from  experimental  data  reported  by  Bjerke 
and  Edmanson  (2004).  All  fractures  were  prohibited  in  simulation  8,  an  illustrative  scenario  designed  for 
comparison  of  results  with  the  other  cases. 

We  prescribed  a  uniform  value  of  Sc=  1.0  pm  for  the  critical  separation  distance  in  the  cohesive  laws  of 
Eqs.  (40)  and  (41),  a  single  value  deemed  best  to  represent  the  effective  fracture  properties  of  the  WHA 
material  system  at  the  length  scale  resolved  by  the  numerical  discretization  (i.e.,  at  the  scale  of  individual 
grains  and  their  interfaces).  Because  these  tungsten  alloys  are  highly  heterogeneous,  with  fracture  behavior 
dictated  by  micro  structural  features  such  as  interfacial  strengths  and  grain  contiguity  and  influenced  by 
large-scale  yielding  (i.e.,  finite  plastic  zones),  difficulties  arise  in  obtaining  consistent  measurements  of 
the  macroscopic  fracture  toughness,  with  values  of  25  <  K[C  <  234  MPa^n  reported  in  the  literature 
(Zamora  et  al.,  1992).  Furthermore,  the  microscopic  fracture  toughness  of  local  interfaces  is  thought  to 
be  smaller  in  magnitude  than  that  of  the  homogenized  material,  as  recovered  uniaxial  test  specimens  reveal 
numerous  microcracks  that  do  not  propagate  to  cause  rupture  (Woodward  and  O’Donnell,  2000).  Assum¬ 
ing  that  the  material  behaves  linear-elastically  and  neglecting  plastic  dissipation  which  may  dominate  en¬ 
ergy  release  in  the  cohesive  zone  (Rice  and  Wang,  1989),  Woodward  and  O’Donnell  (2000)  calculated  a 
static  plane  strain  fracture  toughness  range  of  3.4  ^  Kic  <  7.6  MPav^n  for  microcracking  along  interfaces 
within  a  WHA.  The  mode  I  fracture  energy  corresponding  to  <5c=1.0pm  and  £0  =  2.0  GPa  is 
Gc  =  (l/2)$o<5c  =  1.0  kJ,  which,  upon  assuming  for  illustrative  purposes  only,  an  isotropic  linear-elastic 
stress-strain  response  for  the  composite  micro  structure  with  an  effective  elastic  modulus  of  E—  366  GPa 
and  a  Poisson’s  ratio  of  v  =  0.29  (Dandekar  and  Weisgerber,  1999),  results  in  a  fracture  toughness  value 
of  Klc  =  yjGcE/{\  -  v2)  =  20  MPa-v/5,  thought  here  to  be  a  suitable  compromise  among  the  ranges  re¬ 
ported  in  the  literature.  Note  also  that  our  value  of  5C  is  slightly  larger  than  but  of  the  same  order  of  mag¬ 
nitude  of  the  characteristic  length  of  a  typical  bulk  finite  element  within  the  mesh,  which  we  found  of 
adequate  resolution  to  facilitate  stable,  converged  solutions.  Cohesive  elements  were  deemed  small  enough 
to  resolve  the  process  zone,  the  length  of  which  may  be  liberally  estimated  assuming  a  linear  elastic  bulk 
response  in  the  material  (Rice,  1968;  Espinosa  and  Zavattieri,  2003b).  Specifically,  for  the  present  set  of 
material  parameters,  this  mesh  refinement  criteria  was  satisfied  as  follows: 


7i  E  Gc 

=  2  (T^)^ 


157  pm  »  h  «  0(1  pm), 


(46) 


where  dcz  is  the  cohesive  zone  length  and  h  is  the  typical  element  dimension,  here  on  the  order  of  1  pm. 
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It  should  be  noted  that  the  spall  response  of  the  alloy  is  highly  sensitive  to  heat  treatment.  Dandekar  and 
Weisgerber  (1999)  reported  a  spall  strength  of  1.9  ±  0.4  GPa  for  a  swaged  WHA,  and  a  spall  strength  of 
<1.0  GPa  for  an  annealed  WHA.  The  constitutive  model  parameters  used  for  the  present  set  of  simulations 
correspond  to  the  former  (swaged)  material. 

Fig.  5  shows  the  propagation  and  reflection  of  the  pressure  wave,  and  the  subsequent  spall  fracture,  of 
the  material  aggregate  of  simulation  1.  The  behavior  depicted  in  Fig.  5  is  qualitatively  representative  of  all 
cases  listed  in  Table  3,  except  for  simulation  8,  in  which  fracture  was  prohibited,  and  simulation  9 
(f0  =  2.0  GPa),  in  which  fracture  occurred  prematurely  due  to  resolved  shear  stresses  associated  with  the 
applied  compressive  wave.  Notice  how  the  pressure  pulse  exhibits  modest  heterogeneity  due  to  variable 
properties  of  the  underlying  micro  structure.  By  t  =  35  ns  (Fig.  5(d)),  the  compressive  wave  front  has 
reached  the  free  surface  X1  =  L,  and  tensile  hydrostatic  stresses  (p  <  0)  appear  by  t  =  40  ns,  resulting  in 
spallation  in  the  bulk  of  the  sample  (Fig.  5(f)  and  Fig.  5(g)).  Tensile  pressures  also  arise  on  the  left  side 
of  the  sample  (i.e.,  near  X1  =  0)  in  the  wake  of  the  pressure  pulse,  due  to  the  fixed  boundary  condition 
(43)  applied  for  t  ^  11  ns,  resulting  in  damage  initiation  at  several  locations  towards  this  edge.  Regions 
of  relatively  high  or  low  pressure  do  not  favor  any  one  phase  of  the  microstructure. 

Field  variable  contours  for  simulation  1  are  shown  in  Fig.  6,  prior  to  the  occurrence  of  spall  fracture,  at  a 
simulation  time  of  t  =  20  ns.  The  effective  deviatoric  stress  shown  in  Fig.  6(a)  is  defined  as 

(Te  =  y/3/2(o-pl):(o-pl).  (47) 

We  notice  that  the  distribution  of  effective  stress  ae  is  relatively  heterogeneous  in  the  wake  of  the  applied 
pressure  pulse,  in  contrast  to  the  hydrostatic  pressure  shown  previously  in  Figs.  5(c)  and  (d),  for  example. 
Furthermore,  residual  deviatoric  stresses  remain  throughout  the  deformed  region  of  the  mesh  in  Fig.  6(a),  a 
result  of  incompatible  residual  elastic  deformation  due  to  the  mismatch  in  thermomechanical  properties  in 
different  grains  and  phases  of  the  micro  structure.  Notice  that  deviatoric  stresses  tend  to  concentrate  within 
the  stiffer  W  grains  as  opposed  to  the  more  ductile  matrix  phase.  Effective  plastic  strains  are  shown  in  Fig. 
6(b),  where 

Sp  =  J  v/(2/3)<F  :  <Fd t,  (48) 

with  2<P  =  F  +  Ft  the  rate  of  plastic  stretching  referred  to  the  spatial  frame.  Fig.  6(c)  shows  the  dislocation 
density  pT  introduced  in  Eq.  (25),  a  parameter  which  correlates  strongly  with  the  plastic  strain  of  Fig.  6(b). 
Recall  that  our  variable  pT  represents  the  cumulative  increase  in  shear  strength  at  reference  temperature  60 
averaged  over  all  slip  systems,  and  does  not  account  for  increased  dislocation  mobility  that  would  normally 
accompany  thermal  softening.  From  Figs.  6(b)  and  (c),  respectively,  we  see  that  plastic  strain  and  disloca¬ 
tion  accumulation  concentrate  more  intensely  in  regions  of  the  more  ductile  matrix  phase  of  the  WHA  un¬ 
der  uniaxial  high-strain  rate  compression.  Plastic  strain  magnitudes  within  interior  regions  of  the  mesh 
frequently  exceed  sp  =  0.20,  even  though  the  total  net  strain  of  the  sample  remains  virtually  negligible  at 
the  given  time  instant  of  the  simulation.  Baoping  et  al.  (1994)  measured  a  tendency  for  microhardness  of 
the  binder  phase  to  increase  much  more  intensely  with  increasing  shock  pressure  than  the  microhardness 
of  the  pure  W  grains.  Dislocation  accumulation  was  also  much  more  pronounced  in  the  matrix  phase  rel¬ 
ative  to  the  pure  W  phase  in  the  shock  compression  experiments  of  Baoping  et  al.  (1994),  in  agreement  with 
results  of  our  simulations.  In  Fig.  6(d),  we  see  that  high  temperature  regions  occur  both  in  the  W  grains  and 
in  the  matrix  binder  in  the  wake  of  the  compressive  shockwave.  Notice  also  in  Fig.  6(d)  how  the  temper¬ 
ature  rise  is  greatest  towards  the  left  edge  of  the  sample,  where  the  hydrostatic  stress  state  is  generally  ten¬ 
sile,  and  at  the  lateral  boundaries.  Locally,  temperature  rises  frequently  exceed  several  hundred  Kelvin  in 
our  simulations.  Baoping  et  al.  (1994)  reported  local  temperature  increases  in  excess  of  1000  K  in  a  W  heavy 
alloy,  albeit  shocked  at  much  larger  average  pressures  (in  excess  of  30  GPa)  than  those  imposed  in  the  pres¬ 
ent  simulations.  From  Fig.  6(e),  we  see  that  the  fraction  of  plastic  dissipation  converted  to  heat  energy,  /?  of 
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Fig.  5.  Hydrostatic  pressure  p  (GPa)  for  simulation  1,  at  (a)  t  =  5  ns,  (b)  £=15  ns,  (c)  t  =  25  ns,  (d)  t  =  35  ns,  (e)  £  =  40  ns,  (f) 
£  =  45  ns,  (g)  £  =  50  ns. 


Eqs.  (30)  and  (37),  is  of  generally  lower  value  in  the  pure  W  regions  than  in  the  matrix  regions,  indicating  a 
faster  instantaneous  rate  of  strain  hardening  in  the  former,  at  t  =  20  ns.  However,  values  of  /?  (<0.90)  tend 
to  follow  a  gradient  from  compressed  to  stretched  regions,  with  /?  — ►  1.0  at  the  front  of  the  compressive 
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Fig.  6.  Contours  of  field  variables  for  simulation  1  at  t  =  20  ns:  (a)  effective  stress  oe  (GPa),  (b)  effective  plastic  strain  sp,  (c)  dislocation 
density  pT  (  x  10_7/cm2),  (d)  temperature  6  (K),  (e)  heat  dissipation  parameter  /?,  (f)  recoverable  elastic  energy  density  we  (J/m3),  (g) 
cumulative  plastic  work  wp  (J/m3). 


pulse.  Fig.  6(f)  depicts  the  instantaneous  recoverable  elastic  energy  per  unit  current  configuration  volume, 
i.e., 
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-  \kx&  +  IX?  ■■  ee'|r' ,  (49) 

large  values  of  which  clearly  coincide  with  high  pressure  regions  of  the  matrix  phase,  which  tends  to  exhibit 
relatively  larger  (in  magnitude)  compressive  deformation  d  than  regions  of  the  pure  W  phase,  a  result  of  the 
smaller  effective  bulk  modulus  of  the  former.  The  cumulative  plastic  work  per  unit  current  volume,  defined 
by 

(50) 

is  shown  in  Fig.  6(g),  with  relatively  large  magnitudes  naturally  occurring  in  regions  exhibiting  large  values 
of  temperature  rise  (Fig.  6(d)).  Comparing  Figs.  6(f)  and  (g),  notice  that  the  ratio  of  recoverable  elastic  en¬ 
ergy  to  plastic  work,  w e/wp  1,  in  general.  This  ratio  is  not  to  be  confused  with  the  quantity  1  —  /?,  which 
measures  the  ratio  of  the  residual  lattice  energy  storage  rate  to  the  rate  of  plastic  working  and  generally 
exhibits  values  less  than  unity.  Also,  from  Figs.  6(a)-(g),  we  observe  how  the  rate-  and  history-dependent 
thermoplastic  response  generally  lags  the  instantaneous  elastic  response  (i.e.,  pressure  and  deviatoric  stress 
waves)  at  the  shock  front. 

Average  pressures  p  =  V~l  fpdV  and  temperatures  6  =  V~l  f  OdV  are  shown  in  Figs.  7(a)  and  (b), 
respectively,  for  simulations  labeled  1-8  in  Table  3.  Here  V  is  the  total  reference  volume  of  the  polycrystal¬ 
line  aggregate.  As  noted  previously,  fracture  occurred  prematurely  in  simulation  9,  a  result  of  resolved 
shear  stresses  associated  with  the  applied  compressive  pulse  exceeding  the  prescribed  relatively  low  value 
of  shear  strength  of  cohesive  interfaces,  t0  =  2.0  GPa.  For  simulation  9  the  data  are  not  shown,  as  local 
instances  of  excessive  element  distortion  followed  the  premature  fracture,  necessitating  early  termination 
of  the  simulation  around  t  &  20  ns.  Average  pressure  profiles  in  Fig.  7(a)  are  very  similar  for  all  eight  cases 
prior  to  reflection  of  the  applied  compression  pulse  at  X  =  L,  for  t  <  35  ns.  Peak  pressures  on  the  order  of 
6.0  GPa  occur  at  the  conclusion  of  application  of  the  compressive  pulse  at  t  =  1 1  ns;  subsequently,  the  aver¬ 
age  pressures  decrease  with  time,  in  some  cases  reaching  negative  values  (tension)  upon  wave  reflection  at 
AO  ns,  then  increasing  again  with  time  as  stress  waves  interact  within  the  bulk  material  and  at  the 
boundaries  of  the  domain.  Damage  relaxes  the  reflected  tensile  stresses  in  simulations  1-7,  while  simulation 
eight  supports  greater  tensile  stresses  (peak  of  p  =  —4.62  GPa  at  t  =  66  ns)  since  no  fracture  is  permitted.  In 
simulations  6  and  7,  the  onset  of  fracture  is  delayed  due  to  the  imposed  thermal  strengthening  criterion  at 
interfaces,  thus  explaining  the  larger  peak  average  tensile  stresses  supported  relative  to  simulation  cases  1-5 
for  t  >  50  ns.  Average  temperature  increases  on  the  order  of  400  K  at  t  =  100  ns  are  apparent  in  Fig.  7(b), 
attributed  to  heterogeneous  local  plasticity  and  commensurate  inelastic  dissipation  under  nearly  adiabatic 
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Fig.  7.  Volume-averaged  field  variables  versus  time  t:  (a)  pressure  p,  (b)  temperature  0. 


4632 


J.D.  Clayton  /  International  Journal  of  Solids  and  Structures  42  (2005)  4613-4640 


conditions.  Baoping  et  al.  (1994)  reported  that  temperature  changes  in  excess  of  1000  K  are  expected  in 
tungsten  alloys  subjected  to  shock  pressures  in  excess  of  20  GPa.  Notice  how  the  rate  (slope)  of  temperature 
increase  decreases  slightly  with  time,  as  the  average  stress  state  (and  corresponding  plastic  stress  power) 
relaxes.  Volume-averaged  values  of  the  dissipation  fraction  /?  of  Eqs.  (30)  and  (37)  were  found  to  decrease 
steadily  from  unity  to  approximately  0.70  in  each  simulation,  also  contributing  to  the  slowing  rate  of  tem¬ 
perature  rise  at  later  solution  times. 

Fig.  8  shows  the  total  averaged  free  surface  velocity  histories  x  for  simulations  1-8,  given  by  the  line 
integral 


x 


(51) 


where  the  local  velocity  x1  in  the  integrand  of  Eq.  (51)  is  evaluated  along  the  right  (free)  surface  X1  =  L. 
This  information  corresponds  to  what  one  would  typically  record  in  a  plate  impact  experiment  with  a  veloc¬ 
ity  interferometer  system  (Barker  and  Hollenbach,  1972).  In  general,  comparing  Figs.  7(a)  and  8,  we  see 
that  the  greater  the  intensity  of  damage  within  the  polycrystalline  aggregate,  the  smaller  the  tensile  stress 
and  the  greater  the  average  free  surface  velocity  for  50  ^  f  ^  80  ns,  corresponding  to  the  separation  of  more 
(spalled)  material  from  the  fixed  leftward  portion  of  the  sample.  Maximum  average  velocities  attained 
around  t  =  36  ns  are  of  comparable  magnitudes  among  some  damage  simulations  (cases  5  and  7)  and 
the  damage-free  simulation  (case  8).  For  simulation  8,  average  velocities  oscillate  about  zero  for  t  ^  50  ns, 
as  the  fixed  edge  boundary  condition  at  the  left  side  of  the  mesh,  along  X1  =  0,  prevents  the  polycrystal 
from  undergoing  any  rigid-body  translation.  In  contrast,  for  simulations  1-7,  upon  spallation,  ejected  frag¬ 
ments  of  the  material  are  able  to  continue  to  propagate  away  from  the  bulk  of  the  specimen,  leading  to  a 
significant  positive  free  surface  velocity.  However,  these  ejected  fragments  do  not  regain  the  impact  velocity 
since  the  pieces  are  not  completely  released  in  the  simulations,  as  the  finite  value  of  critical  separation  dis¬ 
tance  specified  for  stress  relief  in  cohesive  zones  (Fig.  2)  is  not  achieved  at  all  interfaces  along  the  spall 
plane.  Differences  in  results  for  simulations  1  and  2  were  negligible  (in  Figs.  7  and  8  as  well  as  in  the  de¬ 
formed  FE  meshes),  as  the  prescribed  shear  strength  of  To  =  4.0  GPa  was  too  large  relative  to  the  tensile 
strength  of  the  interface  (£0  =  2.0  GPa)  to  permit  significant  shear-stress-initiated  fracture.  Comparing  data 
for  simulations  1  and  4,  and  3  and  5,  respectively,  we  notice  effects  of  initial  lattice  orientation.  Differences 
in  predicted  velocities  due  to  micro  structure  are  also  apparent  by  comparing  data  for  simulations  1  (or  2) 
and  3,  corresponding,  respectively,  to  experimental  longitudinal  and  transverse  micro  structural  orientations 
and  grain  shapes  of  Figs.  3  and  4.  Comparing  average  free  surface  velocities  for  simulations  6  and  7  with 
those  for  simulations  1-5,  we  notice  significantly  smaller  velocities  in  the  former  at  later  solution  times 
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Fig.  8.  Spatially-averaged  free  surface  velocity  histories. 
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(t  ^  50  ns).  This  trend  is  due  to  the  relative  reduction  in  spall  separation  associated  with  increased  strength 
with  increasing  temperature  captured  in  simulations  6  and  7.  Upon  consideration  of  Figs.  7  and  8,  we  found 
that  interfacial  properties  (i.e.,  fracture)  apparently  exert  a  more  significant  influence  on  the  averaged  re¬ 
sponse  than  do  bulk  aspects  of  the  micro  structure  (i.e.,  grain  shapes  and  orientations).  Differences  in  micro¬ 
crack  patterns  were  observed  in  deformed  FE  meshes  to  a  degree  commensurate  with  trends  in  the  data  of 
Figs.  7  and  8.  Although  the  properties  of  the  bulk  microstructure  clearly  do  affect  the  free  surface  velocity 
history  (Fig.  8),  their  effects  upon  homogenized  variables  such  as  pressure  and  temperature  (Fig.  7)  are  less 
discernable. 

While  we  could  attempt  to  use  the  data  of  Fig.  8  to  give  some  estimation  of  the  bulk  spall  strength  of  the 
WHA,  the  samples  simulated  contain  too  few  crystals,  and  the  results  too  much  statistical  variability,  to 
make  estimations  that  would  be  worthwhile  to  compare  quantitatively  with  experimentally-obtained  mac¬ 
roscopic  spall  strengths  (Dandekar  and  Weisgerber,  1999).  Furthermore,  the  boundary  conditions  used 
here  are  intentionally  simpler  than  those  typically  encountered  in  flyer-target  plate  impact  experiments, 
to  allow  more  straightforward  interpretation  of  the  numerical  results  with  regards  to  variations  in  material 
properties. 

Fig.  9  gives  the  time  history  of  the  heterogeneity  of  material  (particle)  velocity  along  the  traction-free 
surface  X1  =  L ,  defined  as  the  difference  between  the  local  value  of  i;1  and  the  average  ^  of  Eq.  (51): 


State-of-the-art  methods  are  presently  available  to  achieve  spatial  and  temporal  resolutions  of  particle 
velocity,  in  physical  plate  impact  experiments,  at  length  and  time  scales  comparable  to  those  depicted  in 
Fig.  9  (Vogler  et  al.,  2004).  Fig.  9(a)  gives  data  for  simulation  1,  while  Fig.  9(b)  gives  data  for  simulation 
8.  Recall  that  the  former  simulation  (1)  exhibited  spall  fracture;  the  latter  (8)  did  not.  Grain  shapes  and 
initial  lattice  orientations  were  identical  in  each  of  these  two  simulations.  The  relative  decreases  in  velocities 
near  the  lateral  boundaries  at  X2  =  0  and  X2  =  L  are  attributed  to  attenuation  of  the  pressure  wave  and 
confinement  of  the  mesh  along  these  edges.  For  a  homogeneous  microstructure,  the  free  surface  velocity 
profiles  would  be  symmetric  about  the  centerline  X2  =  L/2  =  75  pm,  and  in  the  absence  of  lateral  boundary 
effects  in  such  a  microstructure,  one  would  have  Ai:1  =  0  for  all  t  ^  0.  Clearly  from  Fig.  9,  neither  profile  is 


Fig.  9.  Free  surface  velocity  perturbation  histories  for  (a)  simulation  1  and  (b)  simulation  8. 
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symmetric  due  to  the  heterogeneous  micro  structure  and  corresponding  variations  in  material  properties 
that  affect  wave  propagation  through  the  polycrystal.  The  profile  in  Fig.  9(a)  is  relatively  more  intense 
and  jagged,  spatially  and  temporally,  than  that  in  Fig.  9(b),  especially  for  t  >  40  ns,  because  of  the  heter¬ 
ogeneous  distribution  of  microcracks  and  fragments  formed  during  spallation  in  later  stages  of  the  former 
simulation. 


5.  Conclusions 

Continuum  crystal  plasticity  models  for  a  two-phase  tungsten  heavy  alloy  have  been  developed,  account¬ 
ing  for  finite  deformation,  plastic  anisotropy,  large  volumetric  strain,  heat  conduction,  thermoelastic  heat¬ 
ing,  rate-  and  temperature-dependent  flow  stress,  and  stored  internal  energy  due  to  lattice  defects.  A  rate 
independent  cohesive  zone  framework  has  been  presented,  capturing  effects  of  temperature  and  mode  mix- 
ity  on  damage  initiation  and  evolution  at  grain  and  phase  boundaries.  Numerical  simulations  of  poly  crys¬ 
talline  aggregates  of  the  alloy  under  shock  compression  conditions  indicate  the  following  trends  in  material 
behavior: 

•  The  pure  W  phase  supports  higher  deviatoric  stresses  than  the  more  ductile  matrix  phase. 

•  The  matrix  phase  accommodates  larger  inelastic  strains  and  exhibits  greater  dislocation  accumulation 
associated  with  strain  hardening,  relative  to  the  W  phase. 

•  The  instantaneous  rate  of  temperature  rise  and  fraction  of  plastic  dissipation  converted  to  heat  energy 
may  favor  either  phase,  depending  upon  the  stress  state  and  time  instant  of  the  deformation  history. 

•  Grain  shapes  and  initial  lattice  orientations  exert  a  noticeable  influence  on  the  shape  and  local  intensity 
of  propagating  and  reflected  stress  waves,  even  in  the  absence  of  fracture. 

•  Interfacial  properties  controlling  the  effects  of  stress  state  (e.g.  mixed  mode  fracture  criteria)  and  temper¬ 
ature  on  microcrack  initiation  and  propagation  significantly  affect  the  evolution  of  spall-induced  dam¬ 
age,  and  hence,  the  overall  spall  strength  of  the  material. 

The  conditions  of  applicability  of  the  model  framework  presented  and  implemented  here  are  now  noted. 
The  thermoelasticity  framework  and  associated  parameters  of  Table  1  are  thought  to  be  reasonably  valid 
for  tungsten  subjected  to  local  pressures  less  than  35  GPa,  corresponding  to  elastic  volume  changes  less 
than  on  the  order  of  10%.  Phase  changes  are  not  captured,  as  transitions  from  BCC  to  FCC  or  DHCP 
(double  hexagonal  close  packed)  structures  are  observed  in  tungsten  at  pressures  on  the  order  of  65  GPa 
(cf.  Ruoff  et  al.,  1998).  Also,  the  melting  temperature  of  pure  W  is  3410  K,  while  that  of  the  matrix  is 
approximately  1750  K  (Zhou  and  Clifton,  1997).  Neither  material  attained  its  melting  point  in  our  simu¬ 
lations.  The  thermoelasticity  model  used  for  the  binder  phase  is  linear  and  isotropic,  due  to  limited  avail¬ 
ability  of  experimental  data.  The  framework  used  for  rate-dependent  plastic  deformation  in  both  phases  of 
the  alloy  (see  also  Clayton,  2005)  was  calibrated  to  the  materials’  homogeneous  stress-strain  responses 
spanning  quasi-static  (0.0001/s)  to  high  strain  rates  (10,000/s),  despite  the  possible  manifestation  of  higher 
local  inelastic  strain  rates  during  simulations  of  wave  propagation  as  occurs  in  plate  impact  experiments. 
The  plasticity  models  may  be  concretely  verified,  and  modified  as  necessary,  for  higher  strain  rate  applica¬ 
tions  (i.e.,  in  excess  of  10,000/s)  as  reliable  experimental  data  become  available.  Additional  avenues  for 
improving  the  physical  description  of  WHA  behavior  in  the  context  of  our  modeling  framework  include 
the  following: 

•  Delineation  of  the  role  of  crystallographic  twinning  observed  in  tungsten  crystals  subjected  to  high  strain 
rates  and  large  pressures  (cf.  Pappu  et  al.,  2001)  from  that  of  dislocation  glide  in  the  plasticity 
framework. 
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•  Characterization  of  non-Schmid  effects  and  slip  strength  anisotropy  on  {112}  planes  as  observed  in 
atomistic  and  experimental  studies  on  pure  W  (Duesberry  and  Vitek,  1998). 

•  Incorporation  of  a  model  for  cleavage  fracture  of  W  grains  along  known  intrinsically  weak  planes  of 
orientation  {100}  and  {110}  and  accounting  for  the  thermal  dependence  of  the  measured  fracture 
toughness  of  these  planes  (Gumbsch  et  ah,  1998). 

•  Depiction  of  length  scale  effects  such  as  grain  size  influences  on  strength,  hardness  and  susceptibility  to 
plastic  flow  instability,  as  experimentally  observed  in  tungsten  during  grain  refinement  (i.e.,  subgrain  for¬ 
mation)  commensurate  with  severe  plastic  deformation  processes  (Aleksandrov  et  al.,  2002). 

We  emphasize  that  the  present  work  offers  an  increased  understanding  of  how  one  may  interpret  exper¬ 
imental  data  on  shock  wave  propagation  and  spall  fracture  in  tungsten  alloys  with  regards  to  effects  of  bulk 
microstructural  and  interfacial  thermomechanical  properties.  Future  studies  will  involve  direct  assessments 
of  results  from  plate  impact  tests  and  numerical  predictions  under  comparable  boundary  conditions.  Spe¬ 
cifically,  comparisons  of  spatially-resolved  particle  velocity  profiles  from  experiments  (Vogler  et  al.,  2004) 
and  simulations  are  envisioned. 
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Appendix  A 


Details  of  the  numerical  procedure  used  for  the  explicit  finite  element  calculations  are  outlined  here.  In 
the  absence  of  body  forces,  the  global  principal  of  virtual  work  corresponding  to  the  linear  momentum  bal¬ 
ance  in  (6)  can  be  written  as  follows  in  the  spatial  frame  (see  e.g.,  Zhou  and  Molinari,  2004): 

f  (px  •  5x)dr  =  —/*(< r  :  5d)dr  —  f  (t  •  [8x])ds,/  +  f  (t  •  5x)dse,  (A.l) 

Jv  Jv  J  Si  J  se 

where  5  denotes  an  admissible  variation,  the  spatial  deformation  rate  tensor  2d  =  |^+  (§^)T,  the  traction 
per  unit  spatial  area  t  =  <x  n,  and  st  and  se  denote  internal  and  external  surfaces  of  the  current  volume  v. 
The  jump  in  virtual  velocity  across  internal  (cohesive)  surface  st  is  written  [5x].  The  energy  balance  (19), 
in  global  variational  form  and  assuming  nonconducting  internal  surfaces,  is  written  as 

yPcmdv= [(tmv- 1  y£d-^)d‘+l  (‘|  ■■««)< >*..  <A-2> 

where 


=  ^(0<)}’(a)  -  p((9{</0  -  0(00{1 k))t  +  p0deeeij/  :  ee. 

a=l 


(A.3) 


Finite  element  discretization  of  (A.l)  leads  to 

MX  =  RV+Ri+Re,  (A.4) 
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where  M  is  the  nodal  mass  matrix,  X  is  the  vector  of  nodal  accelerations,  Rv  is  the  vector  of  internal  forces 
corresponding  to  internal  stresses  and  associated  with  the  first  term  on  the  right  side  of  (A.l),  R1  is  due  to 
forces  exerted  by  cohesive  tractions  associated  with  the  second  term  on  the  right-hand  side  of  (A.l),  and  Re 
is  the  vector  of  externally-applied  forces  resulting  from  third  term  on  the  right-hand  side  of  (A.l).  Analo¬ 
gously,  the  finite  element  equation  corresponding  to  (A. 2)  can  be  written  as  (cf.  Zhou,  1998) 

CE  =  0>  +  KE,  (A.5) 

where  C  is  the  nodal  heat  capacitance  matrix,  E  is  the  vector  of  nodal  temperatures,  includes  the  con¬ 
tribution  of  internal  heating  from  the  first  term  on  the  right  of  (A. 2),  and  K  is  the  matrix  of  heat 
conduction. 

The  solution  algorithm  proceeds  stepwise  in  time  as  follows.  Recall  that  displacements  and  temperatures 
vary  linearly  over  the  3-node  triangular  elements  used  to  model  the  material  continuum.  Let  subscript  n 
denote  the  iteration  number  corresponding  to  the  nth  time  increment  in  the  integration  cycle.  The  acceler¬ 
ation  vector  is  found  from  (A.4)  as 

X„+1  =M-1(R:  +  Rj,  +  R:),  (A.6) 

where  the  right-hand  side  has  been  calculated  from  the  solution  of  the  previous  cycle,  or  from  the  initial 
conditions  for  the  case  of  n  =  0.  The  nodal  velocity  and  position  vectors  are  next  updated  sequentially 
as  follows  (Johnson  et  al.,  1997): 

-^n+l  T  2  (4+1  4—  l)^«+lj  (A. 7) 


XB+1  —  X,  T  (4+ i  4)-^4+i.  (A. 8) 

Analogously,  the  temperature  field  is  updated  explicitly  via 

S„+i  =  Cn  (Ow  +KWE„),  (A. 9) 

^  n+l  ^4  T  (4+1  4)  ^4+1?  (A.  10) 

where  the  temperature-dependent  capacitance  and  conductivity  matrices  have  been  evaluated  during  cycle 
n ,  as  has  the  thermal  force  vector  Q>n,  which  in  our  simulations  consists  primarily  of  the  contribution  to 
heating  from  plastic  dissipation.  Note  also  that  nodal  velocities  and  temperatures  of  (A.7)  and  (A. 9)  are 
updated  consistently  with  externally  imposed  velocity  and  temperature  boundary  conditions  wherever  spec¬ 
ified.  At  this  stage,  local  velocity  gradients  \n+\  are  calculated  at  the  integration  point  (i.e.,  centroid)  of  each 
element,  and  the  deformation  gradient  fw+1  is  found  at  each  integration  point  from  (42),  i.e., 

f«+ i  exp(l„+i(4+i  4))  fn.  (A.  11) 

Correspondingly  the  temperature  is  updated  within  each  element  in  a  manner  consistent  with  Eqs.  (29), 
(31),  (A.9),  and  (A.10): 

0n+ 1  =  On  +  (4+1  —  4  )  @n+\  •  (A.  12) 

The  thermoelastic-viscoplastic  update  within  each  element  then  proceeds  according  to  the  method  dis¬ 
cussed  in  Section  2.5;  for  additional  details  regarding  the  numerical  algorithm,  the  reader  is  referred  to 
McGinty  (2001).  This  constitutive  update,  which  is  the  most  time  consuming  portion  of  the  calculation, 
is  conducted  in  parallel  over  multiple  processors.  The  outcome  of  the  constitutive  update  is  the  local  stress 
field  <rw+i,  which  in  turn  contributes  to  RYn+l  to  appear  in  the  next  iteration  of  (A.6),  as  well  as  the  thermal 
force  contribution  <D„+1  to  appear  in  the  next  iteration  of  (A.9).  We  employ  an  element  failure  criterion  in 
order  to  prevent  excessively  deformed  elements  from  corrupting  the  solution  and  interfering  the  time  incre- 


J.D.  Clayton  /  International  Journal  of  Solids  and  Structures  42  (2005)  4613-4640 


4637 


mentation  scheme  to  be  discussed  shortly  (i.e.,  such  elements  do  not  enter  Eq.  (A.  15)).  Following  common 
practice  in  dynamic  Lagrangian  finite  element  codes  (Johnson  et  ah,  2001),  elements  accumulating  a  very 
large  magnitude  of  plastic  strain,  specifically  £P  >  10.0,  are  not  permitted  to  contribute  to  R^+1,  Ow+1,  and 
Kn+1  for  the  next  cycle,  and  the  local  constitutive  update  is  subsequently  bypassed  in  these  elements.  In  the 
present  context,  this  action  is  justified  by  the  association  of  very  large  inelastic  deformations  with  very  large 
temperature  increases,  resulting  in  relatively  small  stresses  due  to  thermal  softening  (Eq.  (24)).  Masses  and 
capacitances  are  conserved  in  the  nodes  comprising  these  elements,  however.  Recall  also  that  duplicate 
nodes  exist  at  interfaces  between  some  triangular  elements  (in  our  case  along  grain  and  phase  boundaries) 
where  fractures  may  occur.  At  this  point  in  the  computation,  the  normal  traction  £  and  shear  traction  t  (per 
unit  reference  area)  are  computed  along  each  potential  fracture  face  as  an  average  of  the  resolved  stresses 
found  from  the  stress  fields  of  the  two  elements  sharing  that  face.  Criteria  (39)  are  then  tested,  with  the  aver¬ 
age  temperature  at  the  interface  used  when  either  of  s\  or  ti  is  nonzero.  If  neither  initiation  criterion  is 
achieved  at  the  interface,  then  nodal  pairs  at  either  end  that  interface  are  required  to  maintain  the  same 
displacement  and  temperature  values  for  the  next  time  increment,  i.e., 


(Xn+2)\A  =  (Xn+2)\B 


(K+i  +^«+i  +K+i)\a  +  C^h+i  +^«+i  +^»+i)Ib 


Ma  +  Mg 


(A.13) 


(®n+ l)\A 


xi  _  (ffin+l  +^ii+l^»+l)|j  +  (^+1  +^iH-l^»+l)L 
(Ch-OL  +  (Cn+l)\B 


(A.  14) 


for  pinned  coincident  nodes  with  associated  quantities  labeled  by  subscripts  A  and  B.  When  three  or  more 
nodes  are  pinned  at  an  interface  junction  (e.g.,  a  grain  boundary  triple  point),  relations  (A.13)  and  (A.  14) 
are  extended  in  a  straightforward  manner.  On  the  other  hand,  if  either  of  fracture  initiation  criteria  (39)  is 
achieved,  the  nodal  pairs  sharing  the  interface  are  unpinned  and  free  to  acquire  different  positions  and  tem¬ 
peratures  as  the  solution  proceeds.  The  next  term  needed  for  the  subsequent  solution  increment  is  the  force 
vector  R^+1,  which  is  now  computed  from  the  integrating  the  traction  resulting  from  displacement  jumps 
across  active  cohesive  surfaces.  The  traction  contributing  to  R^+1  are  computed  through  relations  (40) 
and  (41),  with  displacement  jumps  calculated  in  terms  of  the  solution  X„+1  from  (A.8).  Additionally,  a 
restoring  force  is  added  to  R^+1  to  prohibit  interpenetration  of  matter  when  compressive  traction  is  trans¬ 
mitted  across  initiated  (unpinned)  cohesive  surfaces.  Finally,  the  contribution  from  external  boundary  con¬ 
ditions,  R®+1,  is  computed.  The  next  integration  cycle  is  then  commenced  (Eqs.  (A. 6)  and  (A. 9)),  with  the 
new  time  increment  chosen  as 


tn+ 2  tn+ 1 


.  (  Qh 

mm  —  .  = 

elements  +  2 [i) j p 


n+l 


(A.15) 


where  q  <  1  is  a  security  coefficient,  h  is  the  element  dimension,  the  first  term  in  parentheses  is  the  constraint 
imposed  such  that  stress  waves  are  properly  resolved  in  the  continuum  elements  (2  is  Lame’s  elastic  stiffness 
constant),  and  the  second  term  in  parentheses  ensures  stable  computations  for  heat  conduction.  The  stress 
wave  constraint  dominated  the  time  incrementation  scheme  for  the  computations  conducted  here  as  a  result 
of  the  specified  mesh  dimensions,  boundary  conditions,  and  material  properties.  When  cohesive  finite  ele¬ 
ments  are  used  in  dynamic  calculations,  a  very  small  value  of  q  is  typically  required  in  order  to  maintain 
stability.  Such  small  time  steps  are  needed  for  two  reasons:  in  order  to  minimize  numerical  smoothing  of 
stress  waves  traversing  a  grid  in  which  duplicate  nodes  are  unpinned  yet  transmit  forces,  and  in  order  to 
ensure  very  small  crack  opening  displacement  perturbations  in  each  iteration  such  that  the  cohesive  trac¬ 
tion-displacement  law  is  adequately  resolved  (Xu  and  Needleman,  1994;  Espinosa  and  Zavattieri, 
2003a, b;  Zhou  and  Molinari,  2004).  In  dynamic  finite  element  simulations  of  cohesive  fracture,  Xu 
and  Needleman  (1994)  and  Zhou  and  Molinari  (2004)  employed  criteria  (A.15)  with  reported  values  of 
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the  security  factor  pranging  from  q  =  0.01  ~  0.1.  In  the  simulations  reported  here,  following  a  trial-and-er- 
ror  investigation,  a  conservative  value  of  q  =  0.01  was  found  to  yield  consistently  converged  results. 
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